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Classical molecular-dynamics simulations were carried out to study epitaxial growth of graphene on 
6H-SiC(0001) substrate. It was found that there exists a threshold annealing temperature above which we 
observe formation of graphitic structure on the substrate. To check the sensitivity of the simulation results, we 
tested two empirical potentials and evaluated their reliability by the calculated characteristics of graphene, its 
carbon-carbon bond-length, pair correlation function, and binding energy. 
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1. Introduction 

Though many experiments were made on various properties of graphene [ij, hmited works have 
been reported on applying the computer simulation technique to study the growth of graphene. One 
very recent work of this kind was done by Tang et al. [2j who used classical molecular dynamics 
(MD) simulation to grow graphene epitaxially on 6H-SiC substrate. The numerical procedure 
which was used for growing graphene on substrate, i.e. by removing layers of Si atoms to mimic 
Si sublimation, is somewhat different from those recent experimental observations by Poon et 
al. [3j who put forth a step-flow growth from edges of terraces. Nevertheless, their simulations 
give a beneficial access to single and multi-layer graphene formation from duly prepared carbon- 
rich layers. Subsequently, simulation studies by the same authors address the thermal stability of 
graphene on 6H-SiC substrate as well [4J, while other simulation works [5, 6j focus more on the 
issue of thermal stability of the buffer carbon layer. As far as classical simulations are concerned, 
the reliability of the simulations depends on the quality of the interaction models, and it is therefore 
of primary importance to assess their capability of treating this system. 

This paper reports a simulation study addressing the following important issues: (a) growing 
graphene epitaxially on the 6H-SiC(0001) substrate, (b) evaluating graphene sheets formed on the 
SiC substrate via comparing two different empirical potentials one of which, the Tersoff poten- 
tial [7], is widely used in the literature and another, its modified version [8j, is supposedly more 
accurate, and (c) discussing the structural stability of mono- and two-layer graphene, based on the 
characteristic features such as the carbon-carbon (C-C) pair correlation function, bond length and 
binding energy. Attempt is made to relate the simulated results to the experimentally observed 
data. 

2. Molecular dynamics method 

2.1. Numerical procedure: LAMMPS software and empirical potentials 

We conduct the MD simulations using the Large-scale Atomic/Molecular Massively Parallel 
Simulator (LAMMPS) software [9j. Two separate MD simulations are carried out, one of which 
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employs the Tersoff potential [7J and the other one uses the modified Tersoff potential [8j (hereafter 
referred to as TEA). It should be noted that while the functional form of both potentials is the 
same, their parameterizations are different, and even their respective range of interaction, that is 
controlled through a smooth cutoff function with two parameters of the potential, are somewhat 
different. We refer the readers to the original papers for more details. Interestingly, the parameter 
files of both these potentials are available in LAMMPS library [9j. To proceed to simulations, we 
first describe how our input data are prepared, and the use of them to study the growth of graphene 
on 6H-SiC(0001) substrate. 



2.2. Numerical procedure: molecular dynamics simulation 



In our Nose-Hoover (NVT ensemble) simulations of the graphene growth, configurations of lay- 
ers of carbon-rich atoms are positioned to loll near a 6H-SiC substrate. Such configurations are 
prepared in the following way. The SiC substrate is set using a crystalline structure with six hexag- 
onal layers repeating periodically in the (0001) direction, each hexagonal layer consisting of two 
sublayers, one for silicon and the other for carbon. Thus, the stacking of 6H-SiC will thus run as AB- 
CACB. . . , which is 6H-SiC(0001). In this work, we focus on the Si-terminated 6H-SiC (see figure [T]) 
and obtain C-rich layers by simply removing the 
topmost Si-layers, which is a numerical proce- 
dure introduced to mimic the sublimation process 
of Si atoms in the epitaxial growth of graphene. 
The dimensions of our orthorhombic simulation 
cell containing the 6H-SiC substrate are 60.07 x 
61.36 X 15.12 in X (generated with lattice pa- 
rameter 2.668 A), y (generated with lattice pa- 
rameter 1.54 A) and z directions, respectively. 
Note that the periodic boundary conditions are 
applied along x and while a vacuum of 30 A IS 
created along the z direction. 

In addition to preparing the initial configu- 
rations of atoms, a technical point which con- 
cerns the growth of multilayer graphene is in or- 
der. In removing silicon atoms directly from the 
6H-SiC crystalline substrate, one focuses on the 
number of C-rich layers by properly choosing pre- 
scribed distances among C-rich layers and sub- 
strate. Consider, for example, four C-rich layers 
after removing Si atoms from the 6H-SiC crys- 
tal. We keep the distance between the substrate 
and the first-layer C-rich atoms next to it at an 
original separation of 1.9 A. Then, between the 
next two, i.e. first and second C-rich layers, we 
make it to lie within 1 A or a separation smaller. 
We set a distance 1.9 A between the second and 
third C-rich layers and resume a separation of 1 A 
again between the third and fourth C-rich layers. 
Note that the C-rich layers take on the original 
crystalline structure with Si atoms removed, i.e. 
a centered hexagonal structure with a C-C bond- 
length of 2.78 A which is longer than 2.65 A of 
Tang et al. [2j. The stringent condition between the C-rich layers and substrate should be strictly 
obeyed otherwise only a few hexagonal rings are created. 

The procedure and the details of parameters used in this simulation growth of graphene on 
6H-SiC substrate are summarized as follows: 
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Figure 1. (Color online) Schematic diagram of 
Si-terminated 6H-SiC(0001). Small blue cir- 
cles are carbon atoms while the large yellow 
ones are Si atoms. 
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• We apply the Tersoff and TEA potentials in two separate simulation series to calculate 
the interatomic C-C and Si-C interactions with parameters given in references ^ and 
respectively. 

• We consider 6H-SiC as a substrate, i.e. six bilayers of Si and C atoms and in each layer of 
Si or C atoms, a total number of 480 atoms is considered. 

• The simulation procedure is performed by first relaxing the system using the conjugate 
gradient minimization method [TU]. The initial distance of C-C atoms in the C-rich layers 
after relaxation is 2.78 A which is larger than the 2.65 A of Tang et al. [2]. Then we perform a 
MD simulation with a timestep At = 0.5 fs, and heat the system until T = 300 K is reached. 
In this process, the temperature is increased using a linear ramp within the Nose-Hoover 
thermostat, so that a heating rate of 10^^ K/s is imposed. At this temperature, the system 
is equilibrated for a time interval of 2 x 10^ At. From this configuration, a set of simulations 
are carried out to increase the temperature of the system to various desired T, in order to 
study the temperature evolution of the carbon layers. As above, a heating rate of 10^^ K/s 
is chosen. At each target temperature T, equilibrium of the system is obtained after a total 
time interval of 3 x 10^ At. In the final stage, the system is annealed down to T = 0.1 K 
at a cooling rate of 5 x 10^^ K/s so that the properties of the carbon sheets can be studied 
without thermal noise. 
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Figure 2. (Color online) Growth of monolayer graphene on 6H-SiC substrate simulated with 
TEA (left) and Tersoff (right) potentials. For TEA potential, the threshold annealing tempera- 
ture at which graphene emerges occurs at a lower temperature within the range 1200 < Ttr < 
1260 K, whereas it is higher within the range 1300 < Ttr < 1450 K for Tersoff potential. 



43802-3 



N. Jakse, R. Arifin, S.K. Lai 



3. Results and discussion 

3.1. Growth of monolayer graphene on 6H-SiC substrate 

The monolayer graphene simulated with the TEA and Tersoff potentials are compared in fig- 
ure |2] The main feature to emphasize is the existence of a threshold annealing temperature Ttr 
signaling the emergence of graphene. For the TEA potential, the carbon layers start to transform at 
1200 K and the graphene layer is formed at Ttr = 1260 K. The latter is close to the experimentally 
value observed by Hannon and Tromp [llj who report the formation of smooth steps of graphene 
in prolonged annealing at 1298 K. For Tersoff 's potential, the transformation of the carbon layers 
to graphene sheet spreads from 1300 K to the threshold value of 1450 K, which is a more gradual 
formation of graphene and at a higher temperature than the experiment [llj. 
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Figure 3. Top row: Comparison of the average C-C bond-length (left) and binding energy 
(right) of monolayer graphene simulated using the Tersoff (open circle) and TEA (solid circle) 
potentials. Bottom row: The pair correlation function of graphene obtained by TEA potential 
(left) at 1260 K and Tersoff potential (right) at 1450 K for a monolayer graphene grown on 
6H-SiC substrate. 



3.2. Bond-length, binding energy, and pair correlation function of monolayer graphene 

The concrete evidence that the TEA potential yields a well-defined graphene structure is its 
prediction of an average C-C bond-length equal to 1.48 A (at 1260 K) (see figures |2]and[3|. This 
value is close to the sp^ -hybridized graphitic carbon (1.42 A). As for Tersoff 's potential, a value of 
1.59 A (at 1450 K) [12j is obtained, which is not as good as the TEA potential. The pair correlation 
functions g{r) in figure [s] are also consistent with the results of bond-length; the position of the 
first maximum of g{r) is 1.487 A for the TEA potential to be compared with 1.508 A for the Tersoff 
potential. Note in figure Is] that the binding energy for the former potential is —7.0941 eV/atom, 
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which is lower than the —6.3620 eV/atom obtained from Tersoff's potential. It is worth mentioning 
that the binding energies have been calculated directly from the potential energy per carbon atom 
given by the Tersoff or TEA potential. Armed with these results, our study for two-layer graphene 
grown on SiC substrate will proceed below by using only the TEA potential. 



3.3. Growth of two-layer graphene on 6H-SiC substrate 

Figure [4] shows the two-layer graphene grown on 6H-SiC substrate. The first graphene layer, 
which is the one that clings to substrate, is relatively less stable than the second layer. This is 
evident from examining Figure [5] for the bond- length and binding energy for the two layers. The 
bond- length of the second graphene layer at T = 1320 K is 1.43 A. This value is smaller than 
1.48 A of the first layer and is closer in magnitude to the 1.42 A bond-length of -hybridized 
graphene carbon. Further evidence can be gleaned also from the pair correlation function g{r) at 
T = 1320 K whose first maximum position for the first layer is 1.53 A to be compared with 1.48 A 
of the second layer. The binding energy of the second layer at T = 1320 K is estimated to be 
—7.1085 eV/atom which is lower than —6.9925 eV/atom for the first layer. 
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Figure 4. (Color online) Two layers of graphene grown on 6H-SiC substrate simulated with TEA 
potential. The first-layer graphene corresponds to the one near the substrate. 
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Figure 5. Average bond-length (right) and binding energy (left) for the first (open circle) and 
second (solid circle) layers graphene grown on 6H-SiC substrate. 



4. Conclusion 

We have applied the MD simulation to the study of the growth of graphene on 6H-SiC substrate 
by epitaxial method. It was found that the choice of an empirical potential in simulation is sensible 
in predicting the threshold temperature at which point the graphene starts to emerge. With the 
TEA potential, we obtained one and two layers of graphene grown on 6H-SiC(0001) substrate. In 
addition to yielding the threshold annealing temperature, which was found to be reasonably close 
to that implied in recent experiments, the characteristics of the grown graphene are confirmed 
by the calculated average C-C bond-length, pair-correlation function and binding energy. With 
this empirical potential in hand, one can envisage investigating the thermal stability of multilayer 
graphene and carrying out a deeper analysis of growth mechanisms. Works along these lines are in 
progress. 
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PicT rpacjDeHy Ha 6H— SiC b ci/iMy/i5ii4i5ix MexoflOM 
MOJieKyji5ipHoY fli/iHaiviiKi/i 

H. >KaKcEl, P. ApiclDeHl^, C.K. /laP 

Hai^ioHa/ibHi/iM no/iirexHinHi/iCi iHCTi/iTyr tpe.Ho6j\9\, VHiBepci/iTeT >Ko3ectDa Oyp'e-Hai^ioHa/ibHi/iCi i^eHTp 
HayKOBi/ix floc/iifl>KeHb, 38402 CaH-MapreH fl'Epe, OpaHij,i5i 

/la6opaTopi5i CK/iaflHi/ix pifli/iH, Oisi/iHHi/iii ctDaKy/ibTer, Hai^ioHa/ibHi/iii i^eHrpa/ibHi/iCi yHiBepci/iTer, Hynr/ii 320, 
TawBaHb 

J\s\9\ Toro, LiJ,o6 flOC/iifli/iTi/i eniraKcia/ibHi/iCi pier rpacfDeHy Ha niflK/iafli^i 6H-SiC(0001) SfliwcHeHO ci/iis/iy/i^iL^iV 
MeroflOM K/iaci/iHHoV MO/ieKy/i5ipHo"i" fli/iHais/iiKi/i. SnaCifleHO icHyBaHH5i nopory reivinepaTypi/i Biflna/iy, Bv\u\e 5ikoV 
cnocTepiraeTbC5i cjDopMyBaHH5i Ha niflK/iaflu,! crpyKrypi/i rpacjDeHy. JA^j\9\ Toro, lu,o6 nepeBipi/iTi/i Hyr/ii/iBicTb 
pesy/ibTariB ci/iis/iy/i^iL^iw, mi/i recTyeis/io flBa eis/inipi/iHHi/ix noreHi^ia/ii/i i oi^inHDeis/io Vxhk) HafliwHicTb lu/i^ixom 
pospaxyHKy xapaKrepi/icri/iK rpacj^eny, flOB>Ki/iHi/i SB'^iSKy Byr/ieUib-Byr/iei^b, napHoV Kope/i^miwHoV cj^yHKi^iV ra 
eHepriV 3B'5i3Ky. 

K/iiOHOBi c/ioea: rpacf^eH, eniraKcia/ibHi/iM pier, MOJieKyjifipHa fli/inaMiKa 
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